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ABSTRACT 


This  tliosis  (lt'\Tlops  ;in  itorativr  algoritfitn  for  the  design  of  AHMA  tiiodcis  of 
signals  in  the  time  floniain.  The  algorithm  is  hasecl  on  optimization  t ('chnicpic's. 
particularly  a  gradient  technique  known  as  the  nslrirtrd  shp  nnthod  is  used.  Phe 
new  algorithm  is  called  the  itcratirf  Hronij  in<  IIkxL  and  the  rc'sults  obtained  using  this 
new  method  are  compared  to  those  obtained  using  the  iterative  preliltering  algorithm. 
The  thesis  shows  that  the  performance  of  the  iterative  Pronv  method  is  in  most  of 
the  cases  comparable  or  superior  to  that  of  the  iterative  preliltering  algorithm. 
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I.  INTRODUCTION 


A.  THE  IDEA  OF  ARMA  MODELING 

The  goal  of  Umar  moth  linfi  to  acniraloly  r('i)n's('iil  an  observed  data  se(|ueiu  (‘ 
as  the  output  of  a  linear  filter.  The  idea  of  representing  a  complicated  procc'ss  with  a 
comparatively  simpler  model  has  many  different  applications.  Curve  fitting  in  mat  h¬ 
ematical  modeling,  analysis  of  electronic  devices  using  eciuivaletit  circuits,  and  system 
transfer  functions  in  automatic  control  are  just  a  few  exampit's  outside  the  area  of  dig¬ 
ital  signal  processing.  I^arametric  modeling  has  also  a  large  number  of  applications  in 
signal  processing.  Currently  there  is  considerable  interest  in  the  parametric  modeling 
approach  to  spectral  estimation.  In  speech  processing  the  applications  include  digital 
transmission,  storage,  and  synthesis  of  the  speech  signal.  Our  particular  interest  is 
the  modeling  of  sonar  signals,  such  as  biologies  and  other  underwater  acoustic  data. 
This  work  forms  part  of  an  overall  research  program  in  sonar  signal  modeling.  The 
research  will  help  to  understand  the  relative  benefits  of  signal  domain  algorithms 
versus  algorithms  based  on  coefficients  of  the  transfer  function.  It  is  hoped  that  the 
method  developed  in  this  thesis  will  become  an  importaiit  tool  in  the  overall  effort 
for  sonar  signal  modeling. 

In  linear  modeling  the  filter  used  to  generate  the  data  secjuence  is  usually  rep¬ 
resented  by  a  linear  difference  equation  with  constant  coefficients.  The  '-transform 
of  this  type  of  system  is  a  rational  polynomial  function,  fhree  type  of  models  arc 
derived  from  this  kind  of  systems;  they  are  known  as  autoregresive  (.\R).  moving 
average  (MA).  and  autoregresive  moving  average  (.XRM.A). 

-Much  work  has  been  done  on  ,\R  models,  which  correspond  to  all-pole  systems. 
The  reason  for  that  is  that  the  parameters  for  the  model  can  be  obtained  by  .solving 


1 


linear  equations,  and  a  great  body  of  theory  has  been  developed  that  applies  to 
this  problem  [Ref.  1.  2].  Relatively  mueh  less  work  has  been  done  with  MA  and 
ARMA  models.  However,  since  M.A  models  have'  limited  applications  and  are  almost 
as  difficult  to  obtain  as  .ARM.-X  models,  most  interest  centers  on  the  latter.  1  he 
hlter  in  this  type  of  models  has  both  poles  and  zc'ros.  and  the  design  lundauK'iitally 
involves  nonlinear  e(|uations.  .\  properly  designed  .ARMA  model  can  provide  better 
performance  than  an  AR  model,  with  a  smaller  number  of  parameters.  .ARMA 
modeling  is  the  topic  of  this  thesis. 

B.  WHY  AN  ITERATIVE  ALGORITHM  IN  THE  SIG¬ 
NAL  DOMAIN 

There  are  many  different  approaches  to  the  problem  of  .ARMA  modeling.  The 
majority  of  them  are  based  on  statistical  techniques  [Ref.  3].  Some  of  these  methods 
regard  the  data  as  a  realization  of  a  random  process  while  others  focus  on  the  data  as 
given  [Ref.  1].  Data  oriented  methods  try  to  minimize  some  criterion  that  estimates 
how  well  the  model  fits  the  data,  in  most  cases  the  least  squares  error  between  the 
signal  and  the  model.  Stochastic  approaches  may  attempt  to  estimate  the  model 
parameters  directly  from  the  data  by  solving  nonlinear  equations  or  by  spectral  fac¬ 
torization.  The  maximum  likelihood  procedure,  for  example  [Ref.  1. 4],  is  essentially 
nonlinear.  A  number  of  indirect  methods  have  been  developed  that  modify  the  norm 
of  the  error  by  separating  the  AR  and  MA  parts  of  the  problem  so  that  at  least  some 
of  the  equations  to  estimate  the  parameters  become  linear.  This  approach  is  found 
in  procedures  such  as  the  Prony's  method.  Shank's  method,  and  the  least  squares 
modified  Yule- Walker  method  [Ref.  1.  2.  5,  6).  .A  different  type  of  approach  replaces 
the  nonlinear  problem  with  iteration  while  trying  to  solve  for  the  .AR  and  MA  pa¬ 
rameters  simultaneously.  The  iterative  prefiltering  method  of  Steiglitz  and  McBride 
[Ref.  7.  8]  is  of  this  type. 
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Our  method  is  also  of  the  latter  type.  However,  the  advantage  is  that  it  works 
directly  with  the  poles  of  the  rational  model,  which  we  know  aliect  the  ptMlortnance 
of  the  system.  The  poles  are  displact'd  in  specific  directiotts  so  that  the  tu'w  model 
minimizes  the  error  between  the  model  output  atid  the  origitial  signal.  Iterative' 
prefiltering  on  the  other  hand  works  with  the  coefficients  of  the  transfer  ftinctioti. 
so  it  is  difficult  or  impossible  to  predict  its  effects  on  the  poles  and  zeros  of  the 
system.  The  new  algorithm  is  much  more  dependable  with  respect  to  convergenct' 
than  the  iterative  prefiltering  algorithm  because  it  moves  [)oles  and  zeros  specifically 
to  minimize  the  error  between  the  model  and  the  original  signal. 

C.  THESIS  OUTLINE 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  II  presents  the 
modeling  methods  that  are  used  in  this  thesis  and  gives  a  brief  explanation  of  all 
of  them.  Chapter  III  introduces  the  reader  to  the  theory  of  multidimensional  op¬ 
timization  by  gradient  methods  and  develops  the  iterative  Prony  method,  which  is 
the  main  contribution  of  this  thesis.  Chapter  IV  presents  the  results  of  testing  the 
algorithm  on  simulated  and  real  acoustic  data  and  compares  these  results  with  those 
obtained  using  iterative  prefiltering.  Chapter  V  gives  conclusions  and  suggestions  for 
future  research. 


II.  ARMA  MODELING  OF  SIGNALS 


A.  MODELING  METHODS  USED  IN  THIS  THESIS 

This  thesis  deals  with  detenninislic  approaciies  to  AKMA  modeling.  1  wo  types 
of  modeling  methods  are  considered.  First  we  have  non-iterative  methods  like  Proriy’s 
method  and  its  alternate  signal  domain  form  [Ref.  1]:  second  we  consider  iterative 
methods  like  iterative  prefiltering  and  the  new  ilcrativt  Pronij  method  developed  in 
this  thesis. 

The  goal  of  Prony's  method  is  to  represent  a  given  sequence  .r[n]  as  the  impulse 
response  of  a  linear  time  invariant  (LTI)  system.  In  the  transform  (c)  domain  this 


representation  has  the  form 


where  X(2)  is  the  ^-transform  of  ifn]  and  8(2)/ A(z)  represents  the  transfer  function 
of  the  system.  This  approximation  is  explained  in  more  detail  in  the  next  section. 
What  has  become  known  as  Prony's  method  in  the  current  signal  processing  literature 
differs  from  Prony's  original  work  in  some  respects.  Our  basic  form  of  Prony's  method 
solves  for  the  coefficients  of  the  transfer  function  in  (2.1). 

The  signal  domain  form  of  Prony's  method  is  closer  to  Prony's  original  work 
[Ref.  1,  p.oGO]  and  seeks  to  represent  the  data  in  terms  of  a  set  of  damped  exponen¬ 


tials  as 


x[n]  s;  nr;*  +  H - h  CpCp 


where  the  are  the  roots  of  the  denominator  polynomial  /!(;)  and  the  Ck  are  the 
complex  coefficients  required  for  the  expansion.  Both  forms  involve  linear  equations 
and  least  squares  techniques. 


■1 


0  ( [/;]  =  .ri//i  —  .r[n] 


Figure  2.1:  Block  diagram  for  tlie  direct  method  lor  signal  modeling. 

The  iterative  prefiltering  method  used  is  due  to  Stciglitz  and  .McBride  [Ref.  7.  S] 
and  attempts  to  match  the  data  with  a  model  of  finite  order  using  an  efficient  iterative 
approach.  It  differs  from  Prony's  method  in  that  it  solves  for  both  numerator  and 
denominator  polynomial  coefficients  ,-iimultaneovsl!j  ixi  each  iteration. 

Whenever  a  signal  is  modeled  using  a  fixed-order  rational  polynomial  model,  an 
approximation  has  to  be  made  and  some  kind  of  measure  has  to  be  used  to  determine 
the  "goodness"  of  the  model.  In  this  thesis  the  least  scpiares  error  norm  {ivboxl^ 
norm)  is  used  to  measure  the  approximation  error.  This  norm  measures  the  energy 
of  the  error  and  is  the  norm  most  widely  used  primarily  due  to  its  mathematical 
tractability  [Ref.  2]. 

B.  OVERVIEW  OF  MODELING  METHODS 
1.  Prony’s  Method 

The  derivations  of  all  the  modeling  methods  presented  in  this  section  follow 
those  in  [Ref.  1.  pp.  .550-564]  and  [Ref.  2].  As  stated  above.  Prony's  method  aims 
at  representing  a  signal  as  the  impulse  response  of  an  LTI  system.  Figure  2.1  shows 
an  implementation  of  this  approximation  which  is  known  as  the  direct  method.  The 
system  function  A(c)  in  the  transform  domain  is  a  rational  polynomial  function 
B{z)lA(z)  with  Q  zeros  and  P  poles.  The  error  signal  ([n]  is  computed  as  the 
difference  between  the  response  of  the  system  .f  [n]  and  the  given  signal  .T[n],  i.e. 

f  [n]  =  .r[n]  —  .f[n].  (‘2.2) 


rile  LTI  system  is  chosen  to  minimize  llie  sum  of  s(|naie(l  (mtois 


(■-'■•5' 

/I. 

This  problem  leads  to  nonlinear  equations  whose  solution  (if  a  unique  solution  actu¬ 
ally  exists)  turns  out  to  be  a  very  difficult  task  [Ref.  1.  2].  do  avoid  these  difficulties, 
a  number  of  indirect  methods  have  been  developed  for  modeling.  Rrony's  method  is 
one  of  these  procedures  and  can  be  derived  as  follows,  hhe  LIT  system  satisfies  the 
difference  equation 

.r[n]  +  aii-[n  -  1]  +  •  •  •  -f  ap.i-[n  -  P\~  -  ■]  + - h  ~ 


If  the  requirement  that 

i[n]  =  j:[n],  n  =  0. 1, - Ns  -  1, 


is  applied  to  (2.4),  where  iV,  is  the  length  of  the  data,  and  the  difference  equation  is 
evaluated  for  n  =  0. 1, . . . ,  .'V’s  —  1,  the  result  is  the  matrix  equation 


■  a-[0] 

0 

0 

...  0 

■  bo  - 

x[l] 

x[0] 

0 

...  0 

b\ 

x[2] 

x[l] 

.r[0j 

...  0 

■  1 

b-i 

a, 

«2 

x[(5] 

x[(5  - 1] 

x[Q  -  2] 

X[Q-P] 

— 

bq 

x\Q  +  l] 

x[(5  -  1] 

■  ■  ■  .c[(5  —  +  1] 

.  . 

0 

x[Ns  —  1] 

■r{Ns  -  2] 

x[Ns  -  3] 

•••  x{Ns-P-[]^ 

.  0 

This  can  be  written  as 


(2.6) 
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where 


X«  = 


,r[0] 

,r[l] 

.r[2l 


0 

■'■[0] 

.r[l] 


0 

0 

.r[0] 


...  u 
0 
u 


and 


X.4  = 


[,r[g]  ,r[g-l]  .v[Q-2\ 
x[Q  + 1]  x[g]  .r[g  -  i] 


AQ  -  A 

.v[Q-P+l] 


(2.7) 


(2.S) 


-1]  .r[iV,-2]  .r[iV,-3]  •••  .r[X..  - -  1] 

and  a  and  b  are  the  vectors  of  coefficients  appearing  in  (2.5).  The  lower  partition 


X  4a  =  0 


(2.!)) 


represents  an  overdetermined  set  of  linear  equations  that  need  not  have  an  e.xact 
solution.  This  set  of  equations  can  be  solved  by  least  squares,  where  a  is  chosen  to 
minimize  the  least  squares  norm  of  the  equation  error  5,4  =  i|e,4j|'^  in 

X,4a  =  e,4.  (2.10) 


This  leads  to  a  set  of  linear  equations  called  the  normal  equations,  which  can  be 
written  compactly  as  [Ref.  1.  pp. 536-537] 


(x*7x,4) 


a  = 


.5,4 

0 


(2.11) 


and  can  be  solved  for  a  and  5,4.  Once  the  vector  a  is  known,  it  can  be  substituted 
in  the  upper  partition  of  (2.6) 

b  =  Xga  (2.12) 


to  solve  for  the  vector  b.  .Although  it  is  referred  here  simply  as  Prony's  method, 
this  procedure  is  also  known  as  the  modern  Prony  method  or  the  e.xtended  Prony 
method. 


I 


x[n 


v(T)  ^-U"]  =  ^  "i"i  -  M”] 


lh)J>\ . /jQ.  0.0.0 - 


Figure  '2.2:  Block  diagram  ior  the  indirect  modeling  problem. 

.Although  Prony's  method  is  simple  to  implement,  it  is  important  to  keep 
in  mind  that  it  is  an  indirect  method.  In  particular  this  procedure  (see  Fig.  2.2) 
minimizes  the  squared  magnitude  of  the  error 


e.4[?7]  =  .r[n]  *  fl[n]  —  b[n] 


(2.1;3) 


where  the  sequences  a[n]  and  6[n]  are  defined  as 


alnl 


On'i  0  <  n  <  P  (Uq  =  1 ) 
0:  otherwise 


and 


,r  1  def  ]  0  <  n  <Q  (fio  =  F 

I  0:  otherwise 


(2.14) 


(2.15) 


where  P  and  Q  are  the  number  of  poles  and  zeros  respectively.  Equation  2.13  repre¬ 
sents  a  different  least  squares  error  from  that  in  (2.3)  where  the  quantity  to  minimize 
was  the  squared  magnitude  of  the  error  e[n]  (see  Fig.  2.1 ).  The  practical  significance 
of  this  difference  is  that  frequently  there  is  a  loss  of  accuracy  in  estimating  the  poles 
and  zeros  by  Prony's  method  [Ref.  3]. 


2.  Signal  Domain  Form  of  Prony’s  Method 


An  alternative  formulation  of  the  method  described  above  can  be  obtained 
if  the  problem  (2.1)  is  stated  in  the  signal  domain  by  representing  .r[n]  in  terms  of  a 
set  of  complex  exponentials; 

x[n]  ~  Cjr"  -|-  +  •  •  •  +  cpCp  (2.16) 
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where  as  stated  before,  the  /‘t  are  the  roots  of  t!ie  polynomial  .l(r).  assumed  to  be 
distinct,  and  the  complex  coelficients  cy.  provi<le  for  t  h('  linear  combination  ol  the  /' 
roots. 

This  ai)proximation  can  be  initially  formulated  as  in  the  s('ction  above  and 
(2.9)  can  be  solved  in  the  least  squares  sense  lor  the  coellicients  ol  .1(c)  (ie.  the 
vector  a).  The  roots  i\.  of  .Kc)  can  then  be  found  and  (2.1b)  can  be  evaluated  for 
n  =  0.  1 . -V,  —  I  to  produce  the  set  of  equations 


This  set  of  equations  can  then  be  solved  in  a  least  stjuares  sense  to  obtain  the  vector 
of  coefficients  c. 

In  the  case  of  multiple  roots  at  the  same  location,  a  slight  variation  of  the 
same  procedure  can  be  used.  Suppose,  for  example,  that  rj  is  a  double  root.  In  this 
case,  the  approximation  is 

x[n]  ^  cir"  +  C2nr'l  H - h  cpr'p  (‘2.18) 

and  the  matrix  equation  to  solve  for  the  coefficients  becomes 


This  situation  is  rare,  however,  because  computational  errors  and  errors  inherent  to 
the  modeling  method  itself  contribute  to  produce  roots  that  may  be  very  close  to 
each  other,  but  not  exactly  at  the  same  location. 
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3.  Iterative  Prefiltering 


The  iterative  |)retilteriri»;  inetluHl  attempts  to  solv(‘  tlie  "din'et"  prohh'iii 
mentioned  iti  snbseetion  1  ol  tins  (  liapt('r  and  to  reiine  t  lie  initial  i)ol<'-z('ro  (  stimate 
by  solving  a  sticcessioti  of  linear  problt'ins.  Idpiatioii  l.'l  (error  for  1  lie  dirc'ct  problem  ) 
can  be  written  in  the  r-domain  as 


E{  =  )  =  X(  =  )  -  X(  =  )  =  .Y(: 


B{z)  X{:)Mz)  -  B{z) 


■  K-] 


-K- 


(2.20) 


The  notion  of  iteration  ran  b('  introduced  that  allows  computation  of  a  new  set  of 
poles  and  zeros  based  on  the  last  known  set  of  poles.  Iterative  prefiltering  replaces 
the  error  for  the  direct  problem  at  tlu'  (/  +  1  iteration  with  the  iterative  error 
function 

If  is  the  inverse  r-transform  of  l//lb)(c).  then  the  error  can  be  written  in  the 

signal  domain  as 


=  j[n]  *  /?b)[n]  *  (2.22) 

The  coefficients  «**'*'**[rt]  and  are  selected  to  minimize  the  corresponding  sum 

of  squared  errors 

=  f2.2;5) 

n  =  H 

at  each  iteration;  this  situation  is  shown  in  Fig.  2.3.  No  general  proof  of  convergence 
has  been  given  for  this  algorithm:  however,  it  is  easy  to  see  that  if  the  iteration  do('s 
converge,  it  must  produce  the  same  answer  as  the  direct  method.  Specifically,  at 
convergence  /fb)  =  and  (2.21)  becomes  the  same  as  (2.20). 

If  we  use  an  indirect  modeling  procedure  like  Prony's  method  to  compute 
the  initial  vector  a  and  we  define  x*') 


xb)|„| 't*'  ^  /dd[i7]. 


(2.21) 
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Fiii,ure  2. IF  Block  (lia!j;raiii  for  t  lie  it('raliv('  prcfiltcriiiii  inclhod. 


then  the  secpieiices  /«'“*[/(]  and  lor  ii  =  0.  1. . .  . ,  A's  -  1  can  l)e  coin[)ule(i  Ironi 

tlie  recursive  dilTerence  e<iuations 


/d‘)[n]  =  /'[n]  -  -  h] 


k=\ 


(2.2d) 


and 


c"M  =  -  *■). 


!.26) 


k=l 


Thus  the  error  (2.22)  can  be  written  as 

p 


-  k]  -  -j]. 


(2.27) 


i:=0 


j=0 


In  order  to  find  the  coefficients  and  the  error  is  written  in  rnatri.x  form 

for  n  =  0.1 . Vj  —  1  as 


X(')  H*’> 


a(«+i) 


=  e‘'+'> 


(2.2S) 


where 


■  X<‘>[P] 

-  1] 

x‘‘)[0] 

X<'>  = 

x<d[p  +  1] 

x“>[li 

.  x(‘)[.'\5-l] 

x'‘)[A^  -  2] 

•••  X<‘)[:V, 

•-!-/>] 

■  /i'‘)[F] 

/!('2[P  _  1] 

...  /d0[p 

= 

1] 

/dd[P] 

-  C?  +  1  ] 

.  1] 

h^%\'s  -  2] 

...  /dO[_v, 

■  -  1  -  f?] 

(2.20) 


(2.20) 


ll 


and 


I'ajuatioii  2. 28  is  aiialuiious  to  (2.1U).  1  Ims  the  least  ■'(jiiari's  iMolileiii  deliiied  1)\ 

iiiiiiiinizintj;  tlu'  norm  i>t'  tlu'  error  in  (2.28)  can  he  redueeil  to 


\vfier<' 

5''""  =  i2.;{:n 

riiis  linear  e((iiati()n  can  tlien  he  .solvt'd  for  (he  xcetors  ol  lilt('r  parameters. 

.\t  this  poitit.  it  is  iri.strucliv<'  to  compare  llu'  p('t formanct'  of  tlie  three 
modeling  nx'thods  outlined  in  this  chai)t(*r  hv  applying  ('ach  to  tnodt'l  a  stnall  segnu'iit 
of  a  transient  sound  corresponding  to  a  wrench  heiiig  struck.  This  wnmch  sound  was 
recorded  and  satnpled  in  the  laboratory  at  a  sampling  rate  of  appro.\imat<'ly  1U.21(J 
Hz.  rids  signal,  denoted  hy  wrenOl.  was  also  used  and  moiieh'd  in  [Ref.  !)].  Fignre 
2.1  shows  a  lUO  point  segmcnit  ot  the  signal  wretiUl  with  one  of  the  thre<'  modc'ls 
(Prony’s,  signal  domain  form  of  Prony's.  and  Iterative  Prelilti'ritig  )  ovi-rlaid.  In  ;dl 
three  <ases  the  signal  was  modeled  with  four  poles  and  four  ziTos.  .\s  can  h(>  st'en. 
the  difference  l)etween  it<‘rative  preliltering  and  tin'  lirst  two  models  is  significant. 
The  non-iterative  methods  match  only  the  initial  points  of  the  se(|uenc('.  and  iirodncc- 
poor  appro.ximations  in  modeling  the  remaining  part  of  the  signal.  I'his  is  due.  in 
large  part,  to  [)oles  not  snfficiently  <lose  to  the  unit  circle.  Iterative  prefilteritig. 
on  the  other  hand,  produces  a  model  which  is  <  lose  to  tlu-  real  se(jU('nc('  alotig  the 
complete  segment. 
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I'igurr  2.4;  Signal  wrc'iiOl  and  its  4  poles/  I  zeros  models,  lal  I  sing  Prony's  method 
(b)  I  sing  Signal  Domain  form  of  Prony’s  metho'l.  |<  i  I  sing  lterati\<'  Pretilterimj,. 


III.  ITERATIVE  ALGORITHM  IN  THE 
SIGNAL  DOMAIN 


A.  MULTIDIMENSIONAL  OPTIMIZATION  BY  GRADI¬ 
ENT  METHODS 

A  multidimensional  function  /(xi .  .rj.  -  ■  ■  .  )  that  is  continuous  and  dilferen- 

tiable  can  be  minimized  using  one  of  s('veral  very  powerful  hillclimh'nuj  techniques 
known  as  gradient  methods  [Ref.  10.  p.  84],  Some  of  those  methods  are  derived  on 
the  basis  of  a  (juadratic  model  that  can  be  obtained  from  a  truncated  Taylor  series 
expansion  of  /(x).  Let  denote  the  value  of  x  at  the  iteration.  Then  for  any 
point  X  =  x'^'*  -f  6  :  when  6  is  small,  the  function  can  be  approximated  by 

=  /<*=>  + (3.1) 

where  g  and  G  represent  the  vector  of  first  derivatives  and  the  matrix  of  second 
derivatives  of  the  function  /(x)  respectively  and  they  should  be  available  at  every 
point.  In  XeiL'ton's  method  the  iterate  is  taken  to  be  x**’*  +  where  the 

correction  minimizes  This  method  is  only  well  defined  when  the  matrix 

of  second  derivatives  G  is  positive  definite,  in  which  case  the  iteration  of  Newton's 
method  is  given  by  the  following  procedure  [Ref.  11.  pp  44-46]: 

1.  solve  =  — g^*"*  for  6  =  A**'* 

2.  set  x*^''"**  =  x**"*  + 

(3.2) 

The  fact  that  G**’’*  may  not  be  positive  definite  when  x*^'*  is  far  from  the  solution, 
and  that  even  when  G*^'  is  positive  definite  convergence  may  not  occur,  makes  this 
method  undesirable  as  a  general  formulation  of  a  minimization  algorithm.  However,  a 


number  of  variations  to  the  basic  method  have  been  proposc'd  that  are  more  suitable' 
for  a  general  class  of  probh'ins.  One  ol  these  methods  is  .\iirtons  milliod  ii'itli  lint 
starch  [Ref.  II.  pp.  17-4f)j  in  which  the  Xewton  algorithm  is  us<’<i  to  generate  a 
direction  of  search 

which  can  later  be  used  in  a  line  search  algorithm  to  actually  calculate  the  correction 
6.  In  the  cases  when  is  not  positive  definite,  the  linear  search  can  be  made 

along  choosing  the  correct  sign  to  ensure  a  descent  direction.  However,  some 

difficulties  that  arise  here  (like  very  high  numerical  costs  and  failure  of  convergence 
for  some  special  cases)  make  this  an  undesirable  approach  for  our  algorithm. 

.As  stated  before.  .Newton's  method  is  defined  only  when  the  matri.v  G***'*  is 
positive  definite,  and  this  matrix  is  positive  definite  only  when  the  error  8  is  "smair’: 
or  better  stated,  the  method  is  defined  only  in  some  neighborhood  Q**'  of  x'*-*  in 
which  agrees  with  +  ^)  in  some  sense.  In  such  cases,  it  is  correct  to 

choose  x*^"^**  =  4-  ^ with  the  correction  ^ minimizing  q^^\6  )  for  all  x<^'^  +  <5 

in  This  method  is  referred  to  as  the  restricted  step  method  because  the  step  is 

restricted  by  the  region  of  validity  of  the  Taylor  series.  [Ref.  11] 

The  region  of  definition  for  the  k'^  iteration  can  be  expressed  as 

=  {x  ;  ||x-x<‘>||  <  (3.4) 

where  ||  ■  ||  denotes  the  norm  of  the  vector.  In  this  case,  the  optimization  problem 
can  be  stated  as: 

minimize^  q^^\8)  subject  to  ||^||  <  (3.5) 

.As  mentioned  before,  the  least  squares  norm  is  the  one  most  commonly  used  in  this 
type  of  problems,  so  it  is  the  one  used  in  this  thesis  and  is  denoted  as  ||  •  II2.  I  he 
problem  that  now  becomes  apparent  is  how  to  select  the  error  margin  of  the 


neighborhood  (d.4).  Fliis  margin  shonid  be  as  large  as  possibh'  Ix'cause  t  lu'  iteration 
step  is  directly  related  to  it.  Various  methods  have  Ixa-n  propost'd  to  control  the 
parameter  one  of  these  methods  attempts  to  insure  that  the  .Newton's  s('arch 

direction  problem  (d.3)  is  always  defined  [Ref.  1 1.  pp.  100-104).  It  does  so  by  adding 
a  multiple  of  the  unit  matri.x  I  to  G**"’  and  computing  the  tiew  problem 


+  //l) 


(4.0) 


where  the  net  effect  is  that  increases  in  cause  ||<5||2  to  decrease,  and  vice  versa. 
If  we  define 


^fik)  f{k)  _  ^ 


(3.7) 


^q(k)  f(k) 

then  the  ratio  represents  a  measure  of  the  accuracy  to  which  approxi¬ 
mates  -f  <5^*’'')  on  the  step,  and  as  the  accuracy  increases  gets  closer  to 


unity.  Using  (3.7),  Marquardt  [Ref.  12]  suggests  an  algorithm  that  tries  to  adaptively 
maintain  as  large  as  possible  while  controlling  the  ratio  The  iteration 
of  such  an  algorithm  is  stated  as: 


1.  given  and  calculate  and  G**’"'; 

2.  factor  G*^'  +  if  not  positive  definite,  reset  i/U)  =  4t/U)  and  repeat: 

4.  solve  (3.6)  to  find  6^^^: 

4.  evaluate  /(x*^*  +  and  hence 


•5.  if  <  0.25  set  =  4i/U4 
else  if  r**'*  >  0.75  set 
else  set 

6.  if  <  0  set  xU'+’>  =  else  set  xU'+')  =  x^d  -|-  ^Ud. 


Here  the  parameters  0.25.  0.75.  4  and  2  are  arbitrary,  and  >  0  is  also  chosen 
arbitrarily  [Ref.  11,  pp.  102-10.4].  Proofs  of  global  and  second  order  convergence  for 
this  algorithm  are  given  in  [Ref.  11.  pp.  96-98]  for  the  cases  when  the  lirst  and  sec¬ 
ond  derivatives  of  the  function  /(x)  exist,  and  the  vector  x^d  belongs  to  a  bounded 
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;)-dimensional  space  for  all  k.  Although  this  method  do(’s  hav('  some  disadvantages, 
it  represents  a  good  basis  for  the  formulation  of  a  general  minimization  algorithm. 
Some  variations  of  this  method  were  rousider<‘d.  but  it  was  tound  that  tor  tin'  s[)e- 
cific  application  of  .\RMA  modeling  in  tlu’  time  domain,  these  variations  [)roduc('<l 
an  extremely  high  overhead  in  calculations. 

B.  THE  ITERATIVE  PRONY  METHOD 

Let  us  now  return  to  the  problem  of  representing  a  sequence'  .r[77]  as  a  linear 
combination  of  complex  exponentials.  Equation  2.17  can  be  written  as 

Rc  =  X  +  c  (d.S) 

where  e  is  called  the  tquation  eiToi\  x  represents  the  data  whicli  may  or  may  not 
be  complex,  c  is  the  vector  of  complex  coefficients,  and  R  is  the  matrix  of  complex 
roots,  which  can  be  written  more  specifically  as 


1 

1 

1 

+ 

j^h 

^R^ 

+ 

Jrh 

•  •  •  rnp 

+ 

jrip 

R  = 

(^Ri 

+ 

jri.f 

{?’R2 

+ 

•••  [t'Rp 

+ 

jripf 

.  ('•Ri 

+ 

{rR2 

+ 

■■■  [r-Rp 

+ 

jripf-'-^  . 

where  r/j^  and  r;,  represent  the  real  and  imaginary  components  of  the  A*’  root,  re¬ 
spectively,  and  Xg  is  the  number  of  data  samples. 

By  defining 

Qllf  ||e||2  =  =  (Rc-xr^(Rc-x).  (3.10) 

it  is  clear  that  the  problem  is  to  find  the  vector  r  =  4-  jr/  of  P  complex  roots  that 

minimizes  the  function  t2(r). 

The  first  and  second  derivatives  of  Q  with  respect  to  r  are  represented  by  the 


vector  g  and  matrix  G. which  are  defined  by 


g  = 


■  iQ 

■  lO 

>'■1, 


:iO 

'"■/p 


(d.ii) 


r 

92  <2 

920 

92  0 

92  0 

:i^Q  1 

OrR^ 

O 

9rR, 

;)2<2 

9rR,  9rRp 

92(2 

,‘^-0 

''•R,  9r/, 

92  0 

,1-^0 

'rRjVrRj 

ilTR^&TRp 

9rR2  9r,| 

■‘TR^arip 

.92  0 

,92  c 

92  0 

92  0 

92  0 

92  0 

drRpdrR^ 

d^O 

dVRpOTR^ 

920 

clTRpdTRp 

a^Q 

9rRp9r,, 
92  0 

aTRpOTJ^ 

9^0 

'5’-RRa>-;p 

92  0 

9’‘/l  9rRj 

92  <1 

9''7,  9rR2 

92(2 

dri^  dTRp 

92  0 

9rf,  9r/j 

92  0 

•)ti  dri^ 

a^<2 

9>'7]  ari 

92  0 

OrI;^^rR^ 

dri^drn^ 

dri^dtRp 

9r/2  9r,j 

ar,^ar,^ 

dri^Orip 

92(2 

92  c 

92  0 

.92  0 

a^Q 

92  0 

.  9r/p9rRj 

9r;p9rRj 

drtpOrRp 

9r/p9r/j 

aripdri^ 

aripdrip  , 

In  order  to  provide  '’or  a  more  compact  notation,  define  the  gradient  operator  with 
respect  to  the  complex  vector  r  as  the  vector  of  partial  derivatives 


r 

2>R, 

9rR2 

- 

9rRp 

Vr,  J 

9 

9rf| 

Orj^ 

L  ^"-'p 

(3.13) 
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consequently.  (3.11)  can  be  expressed  as 


g  = 

Equation  3.12  can  also  be  written  as 


'0 

'O 

>Q 

'>Q 

TO 

'rit 

•"■'l 

■'r/2 

■  1 

■  lO 

'0 

no 

>Q 

no 

TO 

''rn, 

'rf. 

>r,p 

r 

'0 

‘Q 

:iO 

-'0 

1  '"•«) 

■>rr, 

'r,. 

■^rip 

0 

■  rtO 

.TO 

'^0 

:>0 

no 

:dQ 

■a-n. 

■'rpp 

>r,, 

■’rin 

>r,p 

:) 

no 

ao 

no 

no 

no 

.  Otr, 

■^h 

ar,p 

a 

r  .-aT 

;90 

no 

■to 

Or,p 

'•>TRp 

nr,, 

.-Vr;, 

arip 

.\ovv  using  a  somewhat  more  convenient  notation.  (3.15)  can  be  rewritten  as 


G  = 


■) 

f  ^ 

m 

1 

drp^ 

L 

-^r,. 

aQ  ■ 

ar-R, 

nr,^ 

I  =  l...P 
k=\...P 


(3.16) 


and  it  becomes  clear  that  the  matrix  of  second  derivatives  can  be  expressed  compactly 
as 

G=^Vr[(VrQf]=Vr\{Vr,Qf  (Vr^Q)'"'].  (.3.17) 

The  following  subsections  derive  explicit  expressions  for  the  two  quantities  g  and  G. 


1.  Vector  g  of  first  derivatives 


From  (3.8)  it  is  easy  to  see  that 


de  rm 
TT—  =  — c  and 

(Jr,  Or, 


Ori  dr, 
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where  r,  represents  the  real  or  imaginary  parts  of  the  root.  I  sing  (3.18)  and  the 
chain  rule  in  (3.10)  leads  to 


OQ 

f>/?, 


OR 

OrR, 


-c  +  c' 


■dR*^ 


Or 


R, 


(3.10) 


and  explicit  evaluation  of  the  partial  derivatives  of  the  matrix  R  results  in 


OQ 

<^rR, 


0 

1 

+jri.) 

+jriy 


(As  -  1)  (rfi,  ^  J 


cf[0  1  2{rR^-jrjJ  3(rR, -jr/,)'  •••  (A's  -  1)  (fR,  -  jr;,)'^’  ^  ] 


where  c,  is  the  f'*'  component  of  the  vector  c.  If  the  quantity  is  defined  as 


0 

1 


then 


4, 


def 


2(rR.  +  jr  I J 

3(rR,  +jriy 


[N,  -  l)(/x  +jriy 


OR 

Otr^ 


N,-2 


OR- 


Ovr^ 


=  c: 


and  (3.20)  can  be  written  compactly  as 

OQ 


OrR, 


—  _L  c'T 

=  e  ^,  +  c- 


(3.24) 
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At  this  point  it  ran  1)0  n'cognizt'd  that  n'pn'sc-iits  thr  additimi  ol  two  scalar 

([uantities  where  the  lirst  oik*  is  the  complex  coningate  of  the  second,  so  (d.2i)  can 
be  Inrther  simplified  to 

r  7  1 

=  (3.05) 

where  Re[-]  denotes  the  real  part  of  the  vector.  Finally,  using  (d.25)  for  i  —  1  •  ■  ■ 
in  the  upper  partition  of  (d.l  1)  produces 


.\  similar  procedure  can  be  used  to  obtain  the  gradient  of  Q  with  respect 
to  the  imaginary  part  of  the  vector  r.  Once  more,  from  (3. IS)  and  the  chain  rule 
applied  to  (3.10),  it  follows  that 


(3.27) 


These  results  can  be  used  to  generate  an  expression  similar  to  that  in  (3.21)  for  the 
vector  of  first  derivatives  of  Q  with  respect  to  the  imaginary  components  of  r 


(hi^ 


which  in  turn  defines  the  vector  of  partial  derivatives  with  respect  to  the  imaginar\- 
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Equations  .'E26  and  d.dO  can  now  bt'  coinl)ined  as  sliown  in  (d.l  1)  lo  obtain  tlu*  linal 
expression  lor  the  vector  of  (irst  derivativt's 


g  = 


He< 

'  ■ 

r  f  " 

^  1 

\ 

6  > 

-  Cr  . 

' 

Im  < 

C/ 

e  > 

j 

(d.in) 


2.  Matrix  G  of  second  derivatives 


An  expression  for  the  matrix  G  of  second  derivatives  can  be  obtained  as 
follows.  Substituting  (3.24)  and  (3.29)  into  (3.16)  yields 


G 


i  =  I...P 

h=\...P.  (3.32) 


Then  using  the  chain  rule  and  both  expressions  in  (3.18)  leads  to 
G  = 


drn  OrR  Si  ,-,rp  .ir„ 


..7  9C.  ,  ,  ,.-'7nR  ^ 


->rr 


,.7  ''C, 


c-®£i£,-£-'^c-4-. 


•en. 


•"•R, 


■'nu 


:ir,. 


"-n 


I  =  \...P-  k=l...P. 


(3.33) 


From  the  definition  of  ^  in  (3.21 )  it  is  seen  that 


,  _  fief 

Tl -  —  — 


1'-^  ('•/?,  +  ,/'•/.  )■ 


( A.,  -  1 )( .\ .5  -  2)  ( r/},  +  jn,  )  ''■  * 


e,,  if  i=k 


In  the  same  wav  it  can  be  shown  that 


(3.31) 


otherwise. 


-2S.fc 


These  last  four  expressions  together  with  (3.22).  (3.23), (3. 27),  and  (3.28)  can  now  be 
substituted  in  (3.33)  to  obtain 


G  = 


I  =  \...p 


k  =  \...p. 


(3.38) 


Finally  notice  again  that  the  elements  of  the  matrix  G  are  formed  by  addi¬ 
tions  and  subtractions  of  complex  scalars  with  their  respective  complex  conjugates: 
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thus  the  final  expression  for  G  is  given  by 


+  2  Re  [s‘/  cj 

2  Im 

-  2  Im  [s*^^ 

.  Im 

+  2  Im  e] 

2  Re 

-  2  Re  [s^le 

where  it  should  be  remembered  that  s,*..  =  0  for  all  i  7^  k. 

3.  Algorithm  implementation 

An  iterative  method  for  ARM  A  modeling  in  the  time  domain  was  imple¬ 
mented  using  the  results  of  the  last  two  subsections  in  conjunction  with  the  algo¬ 
rithm  presented  in  section  A  of  this  chapter.  We  call  this  method  the  iterative  Prony 
method. 

The  algorithm  uses  the  signal  domain  form  of  Prony ‘s  method  to  calculate 
an  initial  model  for  the  given  sequence.  From  there  it  uses  the  calculated  model 
to  compute  the  error  e.  the  vector  of  first  derivatives  g,  and  the  matrix  of  second 
derivatives  G  and  iterates  until  specific  conditions  are  met.  Figure  3.1  is  an  e.xample 
of  how  the  algorithm  changes  the  position  of  the  poles  and  zeros  of  the  initial  model 
in  order  to  minimize  the  error.  This  figure  represents  the  poles  and  zeros  of  a  transfer 
function  of  order  (-1.3) — 4  poles  3  zeros— that  was  overmodeled  using  a  (6.5)  order 
model.  It  is  clear  from  the  figure  that  the  tendency  in  this  case  is  to  have  a  pole- zero 
cancellation  (see  second  and  third  quadrants)  of  two  poles  and  two  zeros  as  expected. 

Some  features  were  added  to  the  basic  algorithm  in  order  to  deal  with  special 
modeling  cases.  Specifically,  if  the  initial  model  has  some  roots  on  the  real  axis,  then 
because  of  the  way  the  algorithm  iterates,  those  roots  never  move  away  from  the  real 
axis.  A  modification  was  therefore  introduced  to  deliberately  displace  those  roots 
from  the  real  axis  and  proceed  with  the  iterations.  If  the  tendency  of  the  roots  is 
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Figure  3.1:  Displacement  of  the  poles  and  zeros  of  an  iterative  Prony's  model 


“to  go  back”  to  the  real  axis,  then  they  are  returned  to  their  initial  position,  and  the 
iterations  continue.  Otherwise  the  roots  may  continue  to  spread  apart  and  move  as 
a  complex  pair.  Figure  3.2  is  an  example  of  t!ie  displaoe.mcnt  of  the  poles  and  zeros 
of  an  order  (4,3)  model.  In  this  case  the  modeled  signal  actually  has  t,wo  poles  on 
the  real  axis,  and  the  initial  model  correctly  placed  two  of  the  poles  in  the  real  axis. 
Those  poles  are  displaced  from  the  real  axis  by  the  algorithm,  but  then  after  some 
iterations  it  is  clear  that  the  poles  are  tending  to  return  to  the  real  axis.  .At  this  point 
the  poles  are  forced  back  to  the  real  axis  by  setting  their  imaginary  parts  to  zero 
and  the  iterations  continue  until  convergence  is  obtained.  The  opposite  situation  is 
shown  in  Figure  3.3.  In  this  case  the  initial  model  also  has  two  poles  located  on  the 
real  axis,  but  contrary  to  the  case  presented  above,  the  roots,  after  being  displaced 
from  the  real  axis,  continue  to  move  away  from  the  axis  until  they  reach  their  final 
position  in  the  first  and  fourth  cjuadrants  closer  to  the  unit  circle. 
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Figure  3.2:  Displacement  of  the  poles  and  zeros  of  a  1-3  model.  The  modeled  signal 
in  this  case  actually  has  two  poles  on  the  real  a.xis 
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Figure  3.3:  Displacement  of  the  poles  and  zeros  of  a  1-3  model.  Idie  initial  model 
shows  two  poles  in  the  real  axis,  the  final  model  in  this  rase  has  no  poles  in  the  axis 
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IV.  PERFORMANCE  OF  THE  ALGORITHM 
AND  MODELING  RESULTS 


A.  TEST  DATA  USED  IN  THIS  THESIS 

Two  types  of  signals  \v(Me  used  in  tliis  tlu'sis  to  tc'st  tlie  pertortnanee  ol  the 
algorithm.  The  first  type,  whieh  will  he  ealh'd  sminlatid  tc'st  data,  consists  ol  live 
se(|uences  (tOl  to  tOol  each  oiu*  hundred  points  long  that  were  produced  as  the 
impulse  respoi\se  of  a  known  rational  system.  .\ois('  to  produce  an  S.\R  in  the  range 
of  10  to  1")  (IB  was  added  to  the  original  sequence's,  and  the  re’sulting  s('(]uences  were 
designated  as  t01_n  to  t()5_n.  Tlie  origiital  signals  are  d('scrib('d  in  lable  l.l  by  their 
transfer  functions  and  the  location  of  their  poles  and  kteros. 

The  second  group  of  test  signals  consists  of  recorded  nconstic  data.  Two  of  these 
signals  wer(’  recorded  and  sampled  in  the  laboratory.  One  of  them  is  the  sequence 
wrrnOl  already  mentioned  in  (,’hapter  II;  the  other  one  was  obtained  from  human 
speech,  in  particular,  the  signal  voxceLa  corresponds  to  100  samples  of  the  Spanish 
vowel  a.  The  remaining  three  signals  from  the  group  of  acoustic  data  were  recorded 
at  sea  by  a  submarine  platform;  they  correspond  to  sounds  produced  by  marine  life 
and  ice  cracking,  bhe  description  of  the  acoustic  signals  is  ])resent('d  iti  Table  1.2. 

B.  PRESENTATION  OF  RESULTS 

.Ml  simulated  test  signals  were  modeled  twice  with  the  iterative  Prony  tnelhod. 
In  the  first  test  the  e.xact  number  of  poles  and  zeros  of  the  original  model  was  us('d: 
in  the  second  t(’st  all  signals  were  modek'd  using  two  more  pokes  and  zeros  thati  tlu' 
original  model.  This  last  test  is  considered  closer  to  a  real  life  situation  where  the 
exact  order  of  the  signal  to  be  modeled  is  unknown. 
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rABLK  1.1;  nKSClUPTlON  OF  ,s7.\/r/..177;7,)  I  KST  D.VI  A 


POLES 

ZEROS 

().!).')14  Z  X  0.()28(i 

0:  0.770 

().!)  122  Z  ±  0.6200 

0.8087  Z  ±  0.6485 

0;  0.4686 

1.8069:  DC 

0.95  P2  L  ±  1.8846 

0.9514  Z  ±  2.3041 

0;  0.9521 

—0.9521:  DO 

0.9513  Z  ±  0.2097 

0.9049  Z  ±  2.094.3 

0:  DC 

0.9123  Z  ±0.8068 

0.8441;  0.9358 

0;  0.750 

TABLE  4.2:  DESCRIPTION  OF  .4 C’O L'5T/C’ TEST  DATA 


- 1 

SIGNAL 

DESCRIPTION 

wren  01 

Transient  sound  corresponding  to  a  wrench  being  struck 

owe  La 

Spanish  vowel  a 

bio_2133a 

Sperm  whale 

bio_2.385a 

Porpoise  whistle' 

bio_f<0a 

Ice  cracking 

To  mathematically  produce  a  meaningful  measure  of  the  pertormance  ol  a  mod¬ 
eling  algorithm  can  be  quite  difficult  since  various  norms  can  Ix"  deceiving  win'ii 
comparing  errors  of  signals  with  large  differences  in  magnitude.  1  wo  different  ap¬ 
proaches  to  measure  the  performance  of  the  algorithms  were  therefore  used  in  this 
thesis.  The  first  is  quantitative  and  involves  computing  the  squ.ared-norm  of  the  error 
between  the  model  and  the  actual  signal  atid  dividing  it  by  the  total  energy  (norm) 
of  the  signal.  The  second  approach  is  to  overlay  in  a  plot  the  model  and  the  original 
signal  in  order  to  provide  a  visual  comparison  of  the  results.  Phis  is  less  quantitative 
but  frequently  more  revealing  of  errors  in  the  modeling  process. 

1.  Simulated  test  data 

The  first  data  sets  modeled  were  the  simulated  test  data  sets.  Figure  4.1(a) 
is  a  comparison  of  the  normalized  errors  that  result  when  the  sequence  t01_n  is 
modeled  with  2  poles  and  1  zero  using  both  iterative  prefiltering  and  iterative  Proiiy 
methods.  Figure  4.1(b)  and  (c)  show  100  points  of  the  sequence  t01_n  and  the 
two  order  (2,1)  models.  At  this  point  there  is  no  noticeable  difference  between  the 
iterative  prefiltering  and  the  iterative  Prony  models.  Figure  4.2(a)  again  shows  a 
comparison  of  the  normalized  errors  between  an  iterative  prefiltering  model  and  an 
iterative  Prony  model  of  t01_n  for  the  ca.se  wffien  the  signal  t01_n  was  overmodeled 
using  models  of  order  (4,3).  Although  the  difference  between  the  two  modeled  signals 
in  this  case  is  not  large,  notice  that  the  error  for  the  iterative  prefiltering  method 
initially  increases  before  decreasing  while  the  error  for  the  iterative  Prony  method 
decreases  monotonicallv.  This  is  the  first  e.xaniple  of  a  pattern  that  repeats  in  all  but 
one  of  the  simulated  test  signals  that  were  modeled.  Figures  4.3  through  4.10  give 
similar  comparisons  for  the  remaining  simulated  signals.  The  pattern,  which  can  be 
seen  in  Figures  4.1  to  4.10.  is  that  when  the  signals  are  modeled  with  a  number  of  poles 
and  zeros  different  from  that  of  the  actual  order  of  the  signal  (always  overmodeling 
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Figure  4.1:  Signal  tOl-n  and  its  2  poles-1  zero  models,  (a)  .\ornialized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  fOI-ii  and 
the  iterative  prefiltering  model,  (c)  Signal  tOJ-n  and  the  iterative  Prony  model. 
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Figure  4.2;  Signal  tOLn  and  its  4  poles-3  zeros  models,  (a)  .Xormalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  tOl.ri  and 
the  iterative  prefiltering  model,  (c)  Signal  tOl-ii  and  the  iterative  Prony  model. 


.31 


magnitude  magnitude  Error 


10'’ 


10-‘ 


10-2  = 


solid  =  Iterative  Prony  method 
dashed  =  Iterative  prefUtering  method 


5 


_ Iteration  number _ 

10  15  20  25  30 

(a) 


solid  =  signal  t02_n 

5 dashed  =  iterative  prefiltering  model 

A 


0^\ 


\  i 


-5' 


0 


10  20  30  40 


50 

(b) 


60 


70  80  90  100 


solid  =  signal  t02_n 
5  [-  dashed  =  iterative  Prony  model 


,51 - - - - - - - - - , - - - - - - 

0  10  20  30  40  50  60  70  80  90  100 

(c) 


Figure  4.3:  Signal  t02.n  and  its  4  poles-3  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  l02-n  and 
the  iterative  prefiltering  model,  (c)  Signal  t02.n  and  the  iterative  Prony  model. 
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Figure  4.4:  Signal  t02-n  and  its  6  poles-5  zeros  models,  (a)  .Xorrnalized  squared- 
norm  of  the  error  between  the  models  and  the  aetual  signal,  (b)  Signal  t()2^n  and 
the  iterative  prefiltering  model,  (c)  .Signal  t02-n  and  the  iterative  Prony  model. 
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Figure  4.5:  Signal  tOS-Ji  and  its  4  poies-3  zeros  models,  (a)  .Nonnalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  tO-Ln  and 
the  iterative  prefiltering  model,  (c)  Signal  tOi'Ln  and  the  iterative  Prony  model. 
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Figure  4.6:  Signal  tO-Ln  and  its  6  poles-5  zeros  models,  (a)  Normalized  s(|uared- 
norm  of  the  error  between  the  models  and  the  atlual  signal,  (b)  Signal  tOLii  and 
the  iterative  prefiltering  model,  (c)  Signal  iOS-ii  and  the  iterative  Prony  model. 
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Figure  4.7:  Signal  t04-n  and  it.s  4  poles-3  zeros  models,  (a)  .Xormalizcd  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  l()4-n  and 
the  iterative  prefiltering  model,  (c)  Signal  i04-n  and  the  iterative  Prony  model. 
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Figure  4.8:  Signal  t04-n  and  its  6  poles-5  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  l()4-n  and 
the  iterative  prefiltering  model,  (c)  Signal  t04-n  and  the  iterative  Prony  model. 
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Figure  4.9:  Signal  t05-n  and  its  2  poles-1  zero  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t05-n  and 
the  iterative  prefiltering  model,  (c)  Signal  t05-n  and  the  iterative  Prony  model. 
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P'igure  4.10:  Signal  l05-n  and  its  4  poles-3  zeros  models,  (a)  .Xorrnalized  s(|uared- 
norm  of  the  error  between  the  models  and  the  artuai  signal,  (b)  Signal  i()5-n  and 
the  iterative  prefiltering  model,  (r)  Signal  tOo-u  and  the  iterative  Prony  model. 


ill  this  case),  the  behavior  of  tlie  iterative  prelilterine;  method  tc'iids  tcj  deorade  l  i.e.. 
it  *’akes  longer  for  the  algorithm  to  reach  convergenc(')  while  the  behavior  ol  the 
iterative  Prony  method  remains  the  same  or  evi'ii  improves  as  in  the  ease'  of  tOdji 
shown  in  Figures  4.9(a)  and  1.10(a).  Parts  (b)  and  (c)  of  Figure's  1.1  through  4.10 
show  the  original  signals  and  their  respective  models  overlaid.  It  can  be  seen  that 
the  models  arrived  at  by  both  methods  follow  the  original  signals  very  closely  in  all 
cases.  Table  4.3  lists  the  location  of  the  poles  and  zeros  of  the  systems  used  to  model 
all  the  simulated  noisy  signals.  These  systems  were  obtaiiied  using  both  the  iterative 
prefiltering  and  the  iterative  Prony  methods.  The  poles  and  zeros  shown  in  Table  4.3 
can  be  compared  to  the  poles  and  zeros  of  the  original  simulated  signals  presented  in 
Table  4.1.  It  is  clear  that  the  location  of  the  poles  and  zeros  of  the  modeled  signals 
should  not  be  exactly  the  same  as  those  of  the  original  signals  because  some  noise  (in 
the  order  of  10  to  15  dB  SNR)  was  intentionally  added  before  the  modeling  process. 
However,  Tables  4.1  and  4.3  show  a  close  relation  between  the  location  of  the  poles 
and  zeros  of  the  original  sequences  and  the  position  of  the  poles  and  zeros  of  the 
modeled  signals. 

2.  Acoustic  test  data 

As  mentioned  in  the  last  .section  the  acoustic  test  data  represents  sounds 
recorded  both  underwater  and  in  a  laboratory  environment.  In  some  cases  shorter 
segments  were  selected  for  modeling  due  to  the  complexity  of  these  signals.  Once 
again  the  iterative  prefiltering  algorithm  was  used,  and  its  results  were  compared  to 
the  results  obtained  using  the  iterative  Prony  method. 

Before  presentation  of  results,  it  is  important  1o  explain  how  the  model 
produced  by  the  iterative  prefiltering  method  was  selected.  For  all  the  cases,  the 
same  number  of  iterations  was  used  both  for  the  iterative  Prony  and  for  the  iterative 
prefiltering  methods.  In  order  to  obtain  the  best  possible  results  from  the  iterative 
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lABLE  4.3:  POLE-ZERO  LOCATION  OE  HIE  .SVSTE.MS  ESEI)  10  .MODEL 
THE  5/.ULT.1 /’CD  NOISY  TEST  D.VLA 


SIGNAL 
NAME  A- 
(ORDER) 

ITER.YITVE 

PRONY 

ITERATIVE 

PREFILTERING 

POLES 

ZEROS 

POLES 

ZEROS 

t01_n 

(2.1) 

0.9511  L  ±  0.6291 

0:  0.7612 

0.9507  Z  ±  0.6290 

0:  0.7619 

t01_ii 

(4,3) 

0.9510  L  ±  0.6290 
0.9217  L  ±  3.1397 

0.7636 

0.9476;  0.9151 

0.9506  Z  ±  0.6289 
0.9867  Z  ±  2.0944 

0.7630 

0.9886  Z  ^  2.0880 

t02_ii 

(4.3) 

0.9545  Z  ±  0.6238 
0.8455  L  ±  0.6743 

13.2102 

2.1568:  0.1799 

0.9250  Z  ±  0.6313 
0.9142  Z  ±  0.6238 

6.8483 

1.6619:  0.4650 

t02_n 

(6,5) 

0.9541  Z  ±0.6228 
0.8505  Z  ±  0.6797 
0.9896  Z  ±  2.7699 

17.1.527 

2.1237:  0.0952 
0.8760  Z  ±  2.7876 

0.9301  Z  ±  0.6335 
0.9068  Z±  0.6173 
0.9221  Z  ±  2.4.366 

6.5260 

1.8203;  0.5542 
0.9800  Z  ±  2.3896 

t03_n 

(4,3) 

0.9510  Z  ±  1.8850 
0.9505  Z  ±  2.3029 

22.9950 

1.0511;  0.9.365 

0.9499  L  ±  1.8851 
0.9514  Z  ±  2.3030 

12.2468 

0.9680:  0.8020 

t03_n 

(6.5) 

0.9509  Z  ±  1.8851 
0.9506  Z  ±  2.3029 
0.9356  Z  ±  0.2779 

21.1093 

0.9340:  0.6882 
1.1445  Z  ±0.3.507 

0.9492  Z  ±  1.8854 
0.9515  Z  ±2.3028 
0.9796  Z  ±  1.5986 

7.5747 

0.9566:  0.6099 
0.9870  Z  ±  1.6035 

t04_n 

(4.3) 

0.9022  Z  ±  2.0981 
0.9515  Z  ±0.2097 

22.98.53 

0.9494  Z  ±0.8212 

0.9048  Z  ±  2.0956 
0.9516  Z  ±0.2090 

2362.7 

0.9000  Z  ±  0.8090 

t04_n 

(6.5) 

0.9040  Z  ±  2.0966 
0.9517  Z  ±  0.2094 
0.7237:  0.7221 

17..5.349 

0.9104  Z±  0.83.33 
0.7099  Z  ±3.1147 

0.9053  Z  ±  2.0933 
0.9518  Z  ±0.2093 
0.9542;  0.8105 

57.6799 

0.9486  Z  ±  0.8189 
0.9461;  0.7803 

t05_n 

(2.1) 

0.9416:  0.7840 

0.6693 

0.9363;  0.8325 

0.7236 

t05_ii 

(4.3) 

0.9366:  0.8345 
0.9571  Z  ±  2.8523 

0.7330 

0.9571  Z  ±  2.8436 

0.9368:  0.8231 
0.8839  Z  ±  2.7995 

0.7019 

0.8998  Z  ±  2.7977 
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' Ml  tiiD'.i-  I  .iM->  when  ('Diui'i  uciicr  was  luji  obt  aiiird ) . 


'»'.ci’c,i  '::ai  i  •  aTfsiioiiuin'^  !>■<  i  he  Mcraiudi  wiili  the  luwi’st  error 


f  •  ’  f 't  “!  i  '  !  ‘  t *  1(  *1  I !  I  ! 


':,e  Miianiai  ■'e(|ueiiee.  !■  inure  1.11  is  an  exaiii[)le  oi  one  <;1 
■lae  wa'  no!  reaeluai  iisnin  tiu'  iteialivc'  prelilteriiin;  inelhod. 


'oen  ■  !ia! 


'leeted  i>  the  one  obtained  ;il  t  tie  la.st  iteration 


^-'le  .i.neri'iim  luiiudi  in  this  ease  eorresi)t)nds  tt)  a  peak  ot  tlu'  error),  tluni 
■  ne  leNi-.iriTin  iiiociel  oues  !iut  tollows  t  h<' original  sinnal  at  all.  However,  if  we  s('lect 
’.■'.e  -v-reiM  truii;  tiie  iteration  lor  whieh  t he  error  is  the  sniallest  (  17*^‘  iteration  in  this 


M.'e  ,  •  ;ien  v.eoiitain  .i  niodel  that  is  closer  to  the  original  si2;nal.  Some  insight  can 
,;.'o  MO  ootained  '!  \v('  imjk  at  the  Ixdiavior  ol  the  poles  and  zeros  of  the  system  while 
•  :;e  error  ot  tlie  i'er;iti\'e  [tretiltenne  tilu'orithm  is  oscillating,  higures  bl'i  and  4. Id 
'iiow  the  [lositioti  ot  t h<-  [)oh“s  and  zeros  during  t he  oscillation  that  exists  between 
■.•eratii'dis  [.)  fi)  _’(l  ot  the  ca.se  [)resent(>d  in  Figure  1.11.  I'or  the  first  [tart  of  these 
'.•erations  when  rtie  error  i'  low.  the  poles  n'lnain  at  approximately  the  same  position. 
A'  I'eration  _’t).  i  l  in.  l.l  jil  i)  the  poles  suddenly  spread  apart,  some  to  even  outside 
'he  unit  circie,  causing  the  large  increase  in  the<>rror.  .\  similar  t'lfect  occurs  with 
'he  zeros,  although  there  i-  some  significant  movement  of  the  zeros  in  the  t'arlier 
.'eratifdi-..  1  hi.-'  type  ot  behavior  is  avoided  in  the  iterative  Prony  tnethod  where  a 
'ontroileii  and  '’.stematic  displacement  of  the  poles  atid  zeros  is  producc'd  to  finally 
reduce  rhe  error  betweeii  the  mod<'le<l  and  original  signals. 

I  he  sartie  approach  u.>ed  in  the  last  sulrseclioti  was  used  to  model  the  ucoft.s- 
''c  'lata.  .\!l  signals  were  first  mo<leled  with  a  low  order  system  and  subsetjuentlv 
rernofieled  U'lng  a  hiiiher  ordf'r  system.  I  he  results  are  shown  in  h  igures  1. 11  to  4. 23. 
fn  r.'if.isr  of  the  i  a.-^es  the  tiardels  obtained  using  the  it('rative  Prony  method  closelv 
follow  flu  original  'ignah.  ft  can  also  be  observed  that  the  iU'rative  Prony  method 
do#M  fiot  havf  as  large  a  depemlency  in  the  order  ol  the  tnodc'l  seh'ctc'd  as  does  the 
i f era !  1'. e  prch I tfpina  niefhofl.  While  it  is  ol  course  necr-ssjirv  to  seh'ct  a  tnodel  with  a 
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Figure  4.11;  Signal  voweLa  being  modeled  by  iterative  preliltering  using  a  (10.9) 
order  system,  (a)  Normalized  squared-norm  of  the  error  between  the  model  and  the 
actual  signal,  (b)  Signal  vowcLa  and  the  iterative  preliltering  iiiodel  selected  from 
the  20‘^  iteration,  (c)  Signal  voircLa  and  the  iterative  preliltering  model  selected 
from  the  17*^  iteration. 
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Figure  4.12:  Behavior  of  the  poles  of  the  system  used  to  model  the  signal  voweLa 
during  one  oscillation  of  the  iterative  prefiltcring  algorithm,  (a)  15'^  iteration,  (b) 
16*^  iteration,  (c)  17'^  iteration,  (d)  18*^  iteration,  (e)  19"^  iteration,  (f)  20'''* 
iteration. 
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Figure  4.13:  Behavior  of  the  zeros  of  the  system  used  to  model  the  signal  rou'cLa 
during  one  oscillation  of  the  iterative  prefiltering  algorithm,  (a)  15'^  iteration,  (h) 
16*^  iteration,  (c)  17'^  iteration,  (d)  18*^  iteration,  (e)  1!)'^  iteration.  (1)  20'^ 
iteration. 
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high  enough  order  to  properly  model  the  original  signal,  once  such  a  model  has  been 
selected,  increments  or  small  variations  to  its  order  do  not  produce  degradations  on 
the  performance  of  the  algorithm.  On  the  contrary,  it  was  found  that  in  some  cases 
the  performance  of  the  iterative  prefiltering  algorithm  can  be  signilicantly  reduced 
when  the  order  of  the  model  is  increased  slightly.  It  can  also  be  seen  that  the  itera¬ 
tive  Prony  algorithm  tends  to  provide  a  closer  match  to  the  data  than  the  iterative 
prefiltering  method,  and  in  most  of  the  cases  the  rate  of  convergence  of  the  iterative 
Prony  method  was  higher  than  that  of  iterative  prefiltering,  .\nother  important  point 
that  can  be  extracted  from  the  results  presented  in  Figures  4.14  to  4.23  is  that  while 
convergence  with  neither  algorithm  is  guaranteed,  we  obtained  convergence  with  the 
iterative  Prony  method  in  all  cases  for  these  acoustic  signals.  The  same  was  not 
true  for  iterative  prefiltering.  In  most  of  the  cases  the  error  for  the  iterative  Prony 
method  begins  to  decrease  starting  at  the  first  iteration,  and  although  the  change  is 
not  monotonic  in  all  cases,  the  error  after  a  few  iterations  is  consistently  lower  than 
the  initial  error. 
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Figure  4.14:  Signal  wrenOl  and  its  7  poles-6  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  wrenOl  and 
the  iterative  prefiltering  model,  (c)  Signal  wrenOl  and  the  iterative  Frony  model. 
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Figure  4.15:  Signal  wrenOl  and  its  12  poles-11  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
wrenOl  and  the  iterative  prefiltering  model,  (c)  Signal  wrenOl  and  the  iterative 
Prony  model. 
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Figure  4.16:  Signal  roiveLa  and  its  10  poles-9  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  .signal,  (b)  Signal 
voweLa  and  the  iterative  prefiltering  model,  (c)  Signal  rowtLa  and  the  iterative 
Prony  model. 
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Figure  4.17;  Signal  vowcLa  and  its  14  poles-13  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
voweLa  and  the  iterative  prefiltering  model,  (c)  Signal  voweLa  and  the  iterative 
Prony  model. 
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Hgure  4.18:  Signal  bio-2I33a  and  its  8  poles-7  zeros  models,  (a)  .Xormalizcd 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio.2l33a  and  the  iterative  prefiltering  model,  (c)  Signal  bio-2l33a  and  the  iterative 
Prony  model. 
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Figure  4.19:  Signal  hio.2133a  and  its  12  poles-11  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio-2133a  and  the  iterative  prefiltering  model,  (c)  Signal  hio.2133a  and  the  iterative 
Prony  model. 


b3 


magnitude  magnitude  error 


101  - 


IQO  : 


10->  ■ 

10-^  - 
0 


solid  =  iterative  Prony_method 
dashed  =  iterative  prefiltering  method 


5  10  15 

(a) 


Iteration  number 
20  25 


30 


xlO^ 


1  -  solid  =  signal  bio_2385a 
^  dashed  =  iterative  prefiltering  model ' 


Or 

-0.5- 


.11 - - - ^ - 

0  5  10  15  20  25  30  35 

(b) 


n 

40  45  50 


xlO'* 

Ir 

solid  =  signal  bio_2385a 

0.5  ^ 

dashed  =  iterative  Prony  model 

Ob 

> 

-0.5  e'. 

-1L_ 

0  5  10  15  20  25  30  35 

(c) 


r^_ 

40  45  50 


Figure  4.20;  Signal  bio.2385a  and  its  8  poles-7  zeros  models,  (a)  .Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio.2385a  and  the  iterative  prefiltering  model,  (c)  Signal  bio.2385a  and  the  iterative 
Prony  model. 
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Figure  4.21:  Signal  bw.2385a  and  its  12  poles-11  zeros  models,  (a)  .Xormalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (h)  Signal 
bio.2S85a  and  the  iterative  prefiltering  model,  (c)  Signal  bio^2385a  and  the  iterative 
Prony  model. 
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l  igiirc  l.‘2‘2:  Signal  bio.^^On  i\\\(\  its  8  poles-7  zeros  niodi'ls.  (a)  Normalized  s(|uarcd- 
norrn  ol  the  errfir  between  the  models  and  the  actual  signal,  (h)  Signal  bio-SOn  and 
the  iterative  preliltering  model,  (e)  Signal  bio-SOa  and  the  iterative'  Prony  model. 
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Figure  4.23:  Signal  bio.SOa  and  its  12  poIes-11  zeros  models,  (a)  .\ormalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio-SOa  and  the  iterative  prefiltering  model,  (c)  Signal  bio.SOn  and  the  iterative 
Prony  model. 


V.  CONCLUSIONS 


A.  DISCUSSION  OF  RESULTS 

In  this  thesis,  a  new  method  for  modeling  signals  in  the  time  domain  is  developed 
and  applied  to  model  both  recorded  acoustic  data  and  simulated  signals  produced 
as  the  impulse  response  of  a  known  system.  We  call  this  method  the  iterative  Prony 
method.  In  most  of  the  simulated  test  data  sets  the  models  provided  by  the  iterative 
Prony  method  are  sufficiently  close  to  the  original  signals:  in  most  cases,  it  is  diffi¬ 
cult  to  distinguish  between  the  signal  and  the  model.  When  modeling  the  acoustic- 
data  distortion  becomes  apparent  in  some  of  the  models,  which  may  be  due  to  the 
complexity  of  the  structure  of  the  signals.  However,  this  distortion  is  no  worse  than 
for  any  of  the  best  algorithms  that  have  been  used  to  model  this  data  previously. 

The  new  algorithm  was  compared  very  specifically  to  the  iterative  prefiltering 
algorithm  [Ref.  7.  8]  which  has  been  used  in  modeling  a  variety  of  acoustic  data 
[Ref.  3,  9].  The  rate  of  convergence  of  the  iterative  Prony  method  was  in  most  of  the 
cases  comparable  or  superior  to  that  of  the  iterative  prefiltering  algorithm.  Thus, 
while  iterative  prefiltering  sometimes  has  convergence  problems,  the  new  algorithm 
is  much  more  dependable  in  that  respect.  The  price  to  pay  for  this  improvement  is 
in  the  number  of  computations.  While  the  number  of  floating  point  operations  per 
iteration  in  the  iterative  prefiltering  method  is  approximately 

64(P  +  g  -  1  )3  +  SNAP  +  Qf  +  10( P  +  g)A'.  +  12A;, 

iterative  Prony  recjuires  about 

672P^  +  (24W,  +  102)P^  +  (GOA;  -f  46)P  +  198A', 


38 


floating  point  operations  at  oath  ile'ration.  For  example  lor  a  complex  data  st't  ol 
length  100  (A,,  =  100)  and  P  =  Q  —  8  we  have  approximately  lj'2.100  tloatinG: 
point  operations  per  iteration  using  iterative  pretiltering  \ersus  oT'i.OOO  using  iter¬ 
ative  Prony.  If  we  increase  the  order  of  the  model  to  P  =  Q  =  16.  then  we  have 
approximately  2. 787. 824  operations  jter  iteration  in  the  iterative  [)reliltering  algo¬ 
rithm  versus  3.509,560  in  the  iterative  Prony  method. 

B.  RECOMMENDATIONS  FOR  FUTURE  WORK 

The  iterative  prefiltering  algorithm  has  been  the  main  tool  in  the  modeling 
efforts  for  sonar  data  modeling  [Ref.  9.  13].  The  new  iterative  Prony  algorithm  is 
now  at  a  stage  where  it  can  be  substituted  for  the  iterative  prehltering  algorithm  and 
tested  in  operational  use.  To  do  so  needs  some  further  programming  to  make  the 
segmentation  of  the  data  automatic  and  to  make  the  entire  modeling  procedure  more 
of  a  “turn  crank’’  operation.  These  should  be  some  of  the  very  next  steps.  In  addition, 
the  practical  implications  of  the  increased  computation  needs  to  be  addressed,  and 
if  possible  new  methods  need  to  be  developed  to  help  reduce  computations. 

In  a  larger  sense  the  work  reported  in  this  thesis  can  be  used  as  a  base  for 
possible  applications  of  the  iterative  Prony  method  in  the  problems  of  filter  design, 
speech  processing,  and  spectral  estimation.  The  expressions  for  the  vector  of  first 
derivatives  g  and  the  matrix  of  second  derivatives  G  of  the  error  derived  in  Chapter 
III  can  be  used  along  with  different  minimization  methods  to  provide  for  other  new 
modeling  methods  that  may  adapt  better  to  specific  modeling  problems. 
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